Homework 7 Solutions Stat 506 April 2015

  1. Revisit the Oxboys data from HW6.
  1. Refit the random intercept and slope model using a correlation between intercept and slope as in Gelman’s code on p 376. Give appropriate plots and summaries.

No convergence issues, largest R^ is 1.0076644 after 1000 iterations. The jags coefficient and standard deviation estimates are here:

mean sd Rhat n.eff
B.hat[1] 149.3894907 1.7142280 1.006563 380
B.hat[2] 6.5280123 0.3644497 1.005039 500
sigma.b0 8.6347718 1.3826121 1.004578 550
sigma.b1 1.8019277 0.2887957 1.007852 340
rho 0.6077648 0.1406052 1.000374 2000
sigma.y 0.6623925 0.0339128 1.001047 2000

They agree well with lmer estimates:

Std.Dev Correlation
Intercept 8.081 NA
Slope 1.681 0.641
Residual 0.660 NA
Estimate Std. Error t value
(Intercept) 149.371753 1.5854173 94.21605
age 6.525469 0.3363003 19.40370
  1. Do the same using the Inverse Wishart distribution as in his code on p 376–7. Deviances should match, right?

I removed xi.a and got better convergence with max R^ of 1.015. Deviances are pretty close: mean 470.592 (SE = 11.75) for correlation model versus 471.531 (SE 12.169) for the wishart model.

mean sd Rhat n.eff
mu.a 149.3060324 1.5523646 1.000601 5000
mu.b 6.5169483 0.3434120 1.001047 5000
sigma.a 7.9999186 1.1396635 1.000886 5000
sigma.b 1.7210504 0.2592583 1.016492 160
rho 0.6134800 0.1271427 1.010672 260
sigma.y 0.6654763 0.0347933 1.001397 3100
_Again w e have good ag reement in c oefficient and SD estimates._
  1. In the Chapter 15 notes we fit the following model to the owls data.
owls <- read.table("http://www.math.montana.edu/~jimrc/classes/stat506/data/Owls.txt", head=TRUE)
owls$cTime <- owls$ArrvTime - 24
require(lme4)
require(splines)
## Loading required package: splines
owl.fit3 <- glmer(Ncalls ~ offset(log(BroodSize)) + FoodTrt*SexParent +
      SexParent*bs(cTime, df=5) +(1 + cTime|Nest), data=owls, 
      family=poisson)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl =
## control$checkConv, : Model failed to converge with max|grad| = 0.00271455
## (tol = 0.001, component 13)
  1. Translate the model to jags code.
    First I’ll setup the fixed effects part of the model. Log counts of calls are modeled with a log(brrodsize) offset, a treatment, sex-of-parent, treatment by sex interaction, b-spline of degree 5, and interaction between b-spline and sex of parent.
require(R2jags, quietly=TRUE)
cat("model{
  for (i in 1:n){
    y[i] ~ dpois (lambda[i])
    log(lambda[i]) <- offset[i] + X[i,] %*% beta[] + epsilon[i]
    epsilon[i] ~ dnorm (0, tau.epsilon)
  }
  tau.epsilon <- pow(sigma.epsilon, -2)
  sigma.epsilon ~ dunif (0, 100)

  for (col in 1:K){
    beta[col] ~ dnorm (0, .0001)
  }
}
  ",  file = "owlsFE.jags")

logBsize <- log(owls$BroodSize)
X <- with( owls, model.matrix(~ SexParent *(FoodTrt + bs(cTime, df = 5) )))
  ## dim(X)  599 * 14
 owlsData <- with(owls, list("y" = Ncalls, "X" = X, "offset" = logBsize, "n" = nrow(X),"K" = ncol(X) ) )  #, "J" = length(levels(Nest))
  owl.params = c("beta", "sigma.epsilon")  
 owl.jags1 <- jags( data = owlsData, inits = NULL, parameters.to.save = owl.params,  
                 model.file = "owlsFE.jags", n.chains=4, n.iter=1000, progress.bar="none")  
## module glm loaded

Compiling model graph Resolving undeclared variables Allocating nodes Graph Size: 12504

Initializing model

# owl.jags1
 owl.fit1 <- glm(Ncalls ~ offset(log(BroodSize)) + SexParent*(FoodTrt + bs(cTime,5)), data = owls, family=quasipoisson)
  glmCoef <- coef(owl.fit1)
owl.jags1Coef<- apply(owl.jags1$BUGSoutput$sims.matrix, 2, mean)[1:ncol(X)]
owl.jags1SECoef<- apply(owl.jags1$BUGSoutput$sims.matrix, 2, sd)[1:ncol(X)]
  se.lmCoef <- summary(owl.fit1)$coef[,2]
coefTable <- rbind( owl.jags1Coef, glmCoef, owl.jags1SECoef, se.lmCoef)
colnames(coefTable) <- paste("b",1:14,sep="")
knitr::kable(coefTable,digits=3)
b1 b2 b3 b4 b5 b6 b7 b8 b9 b10 b11 b12 b13 b14
owl.jags1Coef 0.724 -1.928 -0.950 0.885 -1.660 -0.081 -1.113 -1.164 0.069 2.618 1.796 3.530 -0.444 3.837
glmCoef 0.618 -1.045 -0.545 1.088 -0.785 0.192 -0.380 -0.565 0.050 1.408 0.827 2.255 -0.830 2.162
owl.jags1SECoef 0.574 1.042 0.174 0.890 0.594 0.916 0.817 0.939 0.230 1.489 0.995 1.456 1.271 1.432
se.lmCoef 0.405 0.729 0.130 0.642 0.422 0.706 0.628 0.727 0.166 1.051 0.688 1.071 0.937 1.011 
owl1.mcmc <- as.mcmc(owl.jags1)
## colnames(owl1.mcmc[[1]])
#for(ndx in 1:length(owl1.mcmc)){
##    owl1.mcmc[[ndx]] <- owl1.mcmc[[ndx]][,c(1, 54, 27, 54:58)]
#}
par(mfrow=c(4,4), mar = c(3,2,2,1))
traceplot(owl1.mcmc)

Having gotten that to run, and seeing fair agreement with a fixed effects glm model, I then add the random nest effects: and intercept, and a slope over centered time.

require(R2jags, quietly=TRUE)
cat("model{
  for (i in 1:n){
    y[i] ~ dpois (lambda[i])
    log(lambda[i]) <- offset[i] + X[i,] %*% beta[] +
         a.raw[nest[i]] + b.raw[nest[i]]*time[i]
    }
  
  for (col in 1:K){
    beta[col] ~ dnorm (0, .0001)
  }
  for(j in 1:J){
    a.raw[j] ~ dnorm(0, tau.a)
    b.raw[j] ~ dnorm(0, tau.b)
    a[j] <- a.raw[j] - mean(a.raw)
    b[j] <- b.raw[j] - mean(b.raw)
  }
  tau.a <- pow(sigma.a, -2)
  sigma.a ~ dunif (0, 100)
  tau.b <- pow(sigma.b, -2)
  sigma.b ~ dunif (0, 100)
  
}
  ",  file = "owlsME.jags")

owlsData2 <- with(owls, list("y" = Ncalls, "X" = X, "offset" = logBsize, "n" = nrow(X),
                            "K" = ncol(X), "J" = length(levels(Nest)), 
                            "nest" = as.numeric(Nest), "time" = cTime ))
owl.param2 = c("beta", "a", "b", "sigma.a","sigma.b")#,"sigma.epsilon")  

 owl.inits <- function (){
  list ( sigma.a = runif(1)*10, sigma.b = runif(1)*10, beta = rnorm(14,0,10))
}

 owl.jags2 <- jags( data = owlsData2, inits = owl.inits, parameters.to.save = owl.param2,  
                 model.file = "owlsME.jags", n.chains=4, n.iter=1000,progress.bar ="none")  
## module glm loaded
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
##    Graph Size: 13806
## 
## Initializing model
 owl.jags2 <- autojags(owl.jags2, n.update = 5, n.iter = 5000,progress.bar ="none")

  1. Run the model to convergence and compare to the glmer output.

The model did converge as maximum R^ is 1.01 and traceplots look good. To see how well lmer and jags agree, I pulled out the fixed effects from each model, and plotted them, then did the same for random nest effects.

The plots above show jags predicted random effects on the horizontal axis against lmer predictions on the y axis. The identity line is shown in grey. One or possibly two nests are a bit above the identity line, but I think this is good agreement.

R/jags Code

Oxboys with correlation.

require(lme4, quietly=TRUE)
require(R2jags, quietly=TRUE)
 data(Oxboys, package="nlme")
 
cat("model {
  for(i in 1:n){
    y[i] ~ dnorm(y.hat[i], tau.y)
    y.hat[i] <- B[Subject[i], 1] + B[Subject[i], 2] * x[i]
  } 
  tau.y <- pow(sigma.y, -2)
  sigma.y ~ dunif (0, 100)
  for (j in 1:J){
    B[j,1:2] ~ dmnorm(B.hat[1:2], Tau.B[,])
  }
  B.hat[1] ~ dnorm(0, .0001)
  B.hat[2] ~ dnorm(0, .0001)
  sigma.b0 ~ dunif (0, 100)
  sigma.b1 ~ dunif (0, 100)
  rho ~ dunif(-1,1)
  Tau.B <- inverse(Sigma.B)
  Sigma.B[1,1] <- sigma.b0^2
  Sigma.B[2,2] <- sigma.b1^2
  Sigma.B[1,2] <- sigma.b0 * sigma.b1 * rho
  Sigma.B[2,1] <- sigma.b0 * sigma.b1 * rho
  
}", 
   file="OB-Corr.jags")

jags2.data = with(Oxboys, list("n" = 234,"J" = 26, "y" = height, "x" = age, "Subject" = Subject)) 
 oxboys.inits <- function (){
  list ( sigma.y = runif(1)*10, sigma.b0 = runif(1)*10, sigma.b1 = runif(1)*4)
}
OB4.params <- c ("B", "B.hat", "sigma.y", "sigma.b0", "sigma.b1", "rho")

OB4.jags <-  jags(data = jags2.data, inits = oxboys.inits, OB4.params,  
                 model.file = "OB-Corr.jags", n.chains=4, n.iter=1000, progress.bar="none")
## OB4.jags
OB4.out <- print(OB4.jags)$summary[c(53:54,57:58, 56, 59), ]
## colnames(OB4.mcmc[[1]])
OB4.mcmc <- as.mcmc(OB4.jags)
 keep <- match(rownames(OB4.out), colnames(OB4.mcmc[[1]]))
for(ndx in 1:length(OB4.mcmc)){
    OB4.mcmc[[ndx]] <- OB4.mcmc[[ndx]][, keep]
}
par(mfrow=c(2,3), mar = c(3,2,2,1))
traceplot(OB4.mcmc) 
maxRhat4 <- max(gelman.diag(as.mcmc(OB4.jags))[[1]][,1])
require(R2jags, quietly=TRUE)
knitr::kable(OB4.out[, c(1:2,8:9)])
require(lme4, quietly=TRUE)  
lmer4 <- lmer(height ~ age + (1 + age|Subject), Oxboys)
## lmer Random Effects:
vc <- VarCorr(lmer4)[[1]]
sigma <- summary(lmer4)$sigma
vc <- matrix(c(attr(vc, "stddev"), sigma,NA,attr(vc, "correlation")[2,1], NA), ncol=2, nrow=3)
dimnames(vc) = list(c("Intercept","Slope","Residual"), c("Std.Dev", "Correlation"))
knitr::kable(vc, digits = 3)

## lmer Fixed Effects:
knitr::kable(summary(lmer4)$coefficients)

Oxboys with Wishart

require(R2jags, quietly=TRUE)
cat("model {
  for(i in 1:n){
    y[i] ~ dnorm(y.hat[i], tau.y)
    y.hat[i] <- a[Subject[i]] + b[Subject[i]] * x[i]
  } 
  tau.y <- pow(sigma.y, -2)
  sigma.y ~ dunif (0, 100)
  for (j in 1:J){
    a[j] <-  B.raw[j,1]
    b[j] <- xi.b* B.raw[j,2]
    B.raw[j,1:2] ~ dmnorm(B.raw.hat[j, ], Tau.B.raw[,])
    B.raw.hat[j, 1] <- mu.a.raw
    B.raw.hat[j, 2] <- mu.b.raw
  }
  mu.a <-  mu.a.raw
  mu.b <- xi.b * mu.b.raw
  mu.a.raw ~ dnorm(0, .0001)
  mu.b.raw ~ dnorm(0, .0001)
  xi.b ~ dunif(0, 100)
  xi.a ~ dunif(0, 100)
  Tau.B.raw ~ dwish(W[ , ], 3)
  Sigma.B.raw <- inverse(Tau.B.raw[,])
  sigma.a <- sqrt(Sigma.B.raw[1,1])
  sigma.b <- xi.b * sqrt(Sigma.B.raw[2,2])
  rho <-  Sigma.B.raw[1,2] / sqrt(Sigma.B.raw[1,1] * Sigma.B.raw[2,2])    
}", 
   file="OB-Wishart.jags")

jags5.data = with(Oxboys, list("n" = 234,"J" = 26, "y" = height, "x" = age, "Subject" = Subject, "W" = diag(2))) 
J <- 26
oxboys.inits <- function (){
  list (B.raw = array(rnorm(2*J),c(J,2)),mu.a.raw = rnorm(1), mu.b.raw=rnorm(1), sigma.y = runif(1)*10, Tau.B.raw = rWishart(1,3,diag(2)))
}
OB5.params <- c ("a","b", "mu.a", "mu.b","rho","sigma.y", "sigma.a", "sigma.b")

OB5.jags <-  jags(data = jags5.data, inits = NULL, parameters.to.save =OB5.params,  
                 model.file = "OB-Wishart.jags", n.chains=4, n.iter=1000, progress.bar="none") 
OB5.jags <-  autojags(OB5.jags,n.iter=5000,Rhat = 1.05, n.update =5, progress.bar="none" )

## OB5.jags

OB5.out <- print(OB5.jags)$summary[c(54,55,57:58, 56, 59), ]
## colnames(OB4.mcmc[[1]])
OB5.mcmc <- as.mcmc(OB5.jags)
 keep <- match(rownames(OB5.out), colnames(OB5.mcmc[[1]]))
for(ndx in 1:length(OB5.mcmc)){
    OB5.mcmc[[ndx]] <- OB5.mcmc[[ndx]][, keep]
}
par(mfrow=c(2,3), mar = c(3,2,2,1))
traceplot(OB5.mcmc) 
mean4 <- round(mean(OB4.jags$BUGSoutput$sims.list[["deviance"]]), 3)
sd4 <- round(sd(OB4.jags$BUGSoutput$sims.list[["deviance"]]), 3)
mean5 <- round(mean(OB5.jags$BUGSoutput$sims.list[["deviance"]]), 3)
sd5 <- round(sd(OB5.jags$BUGSoutput$sims.list[["deviance"]]), 3)
maxRhat5 <- round(max(gelman.diag(as.mcmc(OB5.jags))[[1]][,1]),3)
require(R2jags, quietly=TRUE)
knitr::kable(OB5.out[, c(1:2,8:9)])
owls <- read.table("http://www.math.montana.edu/~jimrc/classes/stat506/data/Owls.txt", head=TRUE)
owls$cTime <- owls$ArrvTime - 24
require(lme4)
require(splines)
owl.fit3 <- glmer(Ncalls ~ offset(log(BroodSize)) + FoodTrt*SexParent +
      SexParent*bs(cTime, df=5) +(1 + cTime|Nest), data=owls, 
      family=poisson)

Owls fixed

require(R2jags, quietly=TRUE)
cat("model{
  for (i in 1:n){
    y[i] ~ dpois (lambda[i])
    log(lambda[i]) <- offset[i] + X[i,] %*% beta[] + epsilon[i]
    epsilon[i] ~ dnorm (0, tau.epsilon)
  }
  tau.epsilon <- pow(sigma.epsilon, -2)
  sigma.epsilon ~ dunif (0, 100)

  for (col in 1:K){
    beta[col] ~ dnorm (0, .0001)
  }
}
  ",  file = "owlsFE.jags")

logBsize <- log(owls$BroodSize)
X <- with( owls, model.matrix(~ SexParent *(FoodTrt + bs(cTime, df = 5) )))
  ## dim(X)  599 * 14
 owlsData <- with(owls, list("y" = Ncalls, "X" = X, "offset" = logBsize, "n" = nrow(X),"K" = ncol(X) ) )  #, "J" = length(levels(Nest))
  owl.params = c("beta", "sigma.epsilon")  
 owl.jags1 <- jags( data = owlsData, inits = NULL, parameters.to.save = owl.params,  
                 model.file = "owlsFE.jags", n.chains=4, n.iter=1000, progress.bar="none")  
# owl.jags1
 owl.fit1 <- glm(Ncalls ~ offset(log(BroodSize)) + SexParent*(FoodTrt + bs(cTime,5)), data = owls, family=quasipoisson)
  glmCoef <- coef(owl.fit1)
owl.jags1Coef<- apply(owl.jags1$BUGSoutput$sims.matrix, 2, mean)[1:ncol(X)]
owl.jags1SECoef<- apply(owl.jags1$BUGSoutput$sims.matrix, 2, sd)[1:ncol(X)]
  se.lmCoef <- summary(owl.fit1)$coef[,2]
coefTable <- rbind( owl.jags1Coef, glmCoef, owl.jags1SECoef, se.lmCoef)
colnames(coefTable) <- paste("b",1:14,sep="")
knitr::kable(coefTable,digits=3)

owl1.mcmc <- as.mcmc(owl.jags1)
## colnames(owl1.mcmc[[1]])
#for(ndx in 1:length(owl1.mcmc)){
##    owl1.mcmc[[ndx]] <- owl1.mcmc[[ndx]][,c(1, 54, 27, 54:58)]
#}
par(mfrow=c(4,4), mar = c(3,2,2,1))
traceplot(owl1.mcmc)

Owls Mixed

require(R2jags, quietly=TRUE)
cat("model{
  for (i in 1:n){
    y[i] ~ dpois (lambda[i])
    log(lambda[i]) <- offset[i] + X[i,] %*% beta[] +
         a.raw[nest[i]] + b.raw[nest[i]]*time[i]
    }
  
  for (col in 1:K){
    beta[col] ~ dnorm (0, .0001)
  }
  for(j in 1:J){
    a.raw[j] ~ dnorm(0, tau.a)
    b.raw[j] ~ dnorm(0, tau.b)
    a[j] <- a.raw[j] - mean(a.raw)
    b[j] <- b.raw[j] - mean(b.raw)
  }
  tau.a <- pow(sigma.a, -2)
  sigma.a ~ dunif (0, 100)
  tau.b <- pow(sigma.b, -2)
  sigma.b ~ dunif (0, 100)
  
}
  ",  file = "owlsME.jags")

owlsData2 <- with(owls, list("y" = Ncalls, "X" = X, "offset" = logBsize, "n" = nrow(X),
                            "K" = ncol(X), "J" = length(levels(Nest)), 
                            "nest" = as.numeric(Nest), "time" = cTime ))
owl.param2 = c("beta", "a", "b", "sigma.a","sigma.b")#,"sigma.epsilon")  

 owl.inits <- function (){
  list ( sigma.a = runif(1)*10, sigma.b = runif(1)*10, beta = rnorm(14,0,10))
}

 owl.jags2 <- jags( data = owlsData2, inits = owl.inits, parameters.to.save = owl.param2,  
                 model.file = "owlsME.jags", n.chains=4, n.iter=1000,progress.bar ="none")  
 owl.jags2 <- autojags(owl.jags2, n.update = 5, n.iter = 5000,progress.bar ="none")
owl2.mcmc <- as.mcmc(owl.jags2)
## colnames(owl2.mcmc[[1]])
for(ndx in 1:4){
    owl2.mcmc[[ndx]] <- owl2.mcmc[[ndx]][,c(1, 27:28, 54:55, 68:71  )]
}
par(mfrow=c(3,3), mar = c(3,2,2,1))

traceplot(owl2.mcmc)
maxOwlRhat  <- round(max(gelman.diag(as.mcmc(owl.jags2))[[1]][,1]),3)

Owl Plots

require(lme4, quietly=TRUE)
require(R2jags, quietly=TRUE)
par(mfrow=c(1,3), mar = c(2,2,2,1))
owl.jags2Coef<- apply(owl.jags2$BUGSoutput$sims.matrix[,55:68], 2, mean)[1:ncol(X)]
owl.jags2SECoef<- apply(owl.jags2$BUGSoutput$sims.matrix[,55:68], 2, sd)[1:ncol(X)]
  glmerCoef <- summary(owl.fit3)$coef[,1]  
  se.glmerCoef <- summary(owl.fit3)$coef[,2]
plot(owl.jags2Coef, glmerCoef, main = "Fixed Effects")
abline(0,1, col = "grey")
 glmerREs <-  ranef(owl.fit3)[[1]]
 jagREs <- summary(as.mcmc(owl.jags2))$statistics[1:54,1]
plot(jagREs[match( paste("a[",1:27,"]",sep=""),names(jagREs)[1:27])], glmerREs[,1], main = "Random Intercepts")
abline(0,1, col = "grey")
 plot(jagREs[match( paste("b[",1:27,"]",sep=""),names(jagREs)[1:27+27])+27], glmerREs[,2], main = "Random Slopes")
abline(0,1, col = "grey")